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Abstract 

In this paper we present a new approach towards global passive approximation in order to find a passive 
transfer function G(s) that is nearest in some well-defined matrix norm sense to a non-passive transfer 
function H(s). It is based on existing solutions to pertinent matrix nearness problems. It is shown that the 
key point in constructing the nearest passive transfer function, is to find a good rational approximation of 
the well-known ramp function over an interval defined by the minimum and maximum dissipation of H(s). 
The proposed algorithms rely on the stable anti-stable projection of a given transfer function. Pertinent 
examples are given to show the scope and accuracy of the proposed algorithms. 
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1. INTRODUCTION 

For linear time-invariant systems, passivity guarantees stability and the possibility of synthesis of a 
transfer function by means of a lossy physical network of resistors, capacitors, inductors and transformers 
[l[ . Therefore, passivity enforcement Q and passification (passivation) Q have become important issues in 
recent years (j48| > especially as more and more software tools render transfer functions which need passivity 
enforcement as a postprocessing step in order to generate reliable physical models. However, most of the 
techniques 0-0] are local perturbative and/or feedback approaches with fixed poles, while [1] is based on 
Fourier approximation, yielding passivated systems with a large number of poles. 

In this paper we present a new global approach in the sense that we find a passive transfer function G(s) 
that is nearest in a well-defined matrix norm sense to a non-passive transfer function H(s). It is based 
on existing solutions to some pertinent matrix nearness problems 0, [Io| . We show that the key point in 
constructing the nearest passive transfer function G(s), is to find a good rational approximation for the ramp 
function max(0, x) over an interval defined by the minimum and maximum dissipation of the non-passive 
transfer function H(s). It is also shown that in the Chebyshev or minimax sense this requires finding a 
rational Chebyshev approximation of the square ro ot ~J x over the interval [0,1]. The proposed algorithms 
rely heavily on the stable anti-stable projection [ll|, [l2| of a given transfer function. Finally, five pertinent 
examples, both SISO and MIMO, are given to show the accuracy and relevance of the proposed algorithms. 



2. PASSIVITY AND DISSIPATION 

Notation : Throughout the paper X T and X H respectively denote the transpose and Hermitian transpose 
of a matrix X, and /„ denotes the identity matrix of dimension n. The Frobenius norm is defined as 
= VtiX^X and the spectral norm (or 2-norm or maximum singular value) is defined as ||X||2 = 
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\J ^m&x{X H X). It is easy to show that H-X^H^ = ||^||f and ||X ff ||2 = ||-X"||2- For two Hermitian matrices 
X and Y, the matrix inequalities X > Y or X > Y mean that X — Y is respectively positive definite or 
positive semidefinite. The closed right halfplane 5fte [s] > is denoted C+. 
For the real system with minimal realization 



= Ax 
= Cx 



Bu 
Du 



(la) 
(lb) 



where B ^ 0, C ^ are respectively n x p and p x n real matrices and A ^ is a n x n real matrix, to be 
passive, it is required that the p x p transfer function 

H(s) = C(sl n — A)~ 1 B + D 

is analytic in C+, such that 

H(iu) + H(ico) H > VweK 

It is well-known that the positive-real lemma in linear matrix inequalty (LMI) format : 3 P T = P > 
such that 



A T P + PA 
B T P-C 



PB 
—D 



C T 
D T 



< 



guarantees the passivity of the system ([1}. A necessary, but not sufficient, condition for passivity is that A 
is stable, i.e., its eigenvalues are located in the closed left halfplane. In the sequel we will always suppose 
that A is Hurwitz stable, i.e., its eigenvalues are located in the open left halfplane. We will also assume, 
unless otherwise stated, that H(s) is non-passive, and devise ways of finding another as close as possible 
passive transfer function G(s). 

In order to measure how far a given system is from passive we define the minimum dissipation 8- (H) [Tij ] 
as 



8-(H)=udnX^[R(u)] 



where 



R(u) = H(iuj) + H(iuj) H 
Similarly, we also define the maximum dissipation 6+(H) as 

5+(H) = maxA max [R(lo)} 

It is clear that the system is passive if and only if 5 -(H) > 0. If S-(H) < the system is non-passive, and 
if 8+(H) < 0, the system is anti-passive, in the sense that then the system with transfer function —H(s) is 
passive. 

In the sequel we will assume, unless otherwise stated, that the system is non-passive but passifiable, i.e., 
— oo < S -(H) < < 5 + (H) < oo. To obtain 5 -(H) (or similarly 5+(H)), a simple bisection algorithm, based 
on the existence (or non-existence) of imaginary eigenvalues of the one-parameter Hamiltonian matrix 



" A 







B 





—A T _ 


+ 





(5Jp -D- D?)- 1 [C B T ] 



was proposed in We have 

Proposition 2.1. S > S-(H) if and only ifN$ admits purely imaginary eigenvalues. 



Proof. See [14 1 . 



□ 



It is clear that Proposition 12.11 always allows to decide, by checking the eigenvalues of IM5, whether 
8 > 8 -(H) or not. This forms the basis of the bisection algorithm of [13]. The only problem is to start with 
a so-called bracket, i.e., provable lower and upper bounds for 8-(H). For that purpose we have 
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Proposition 2.2. 

- 2||#|| 00 < S-(H) < A min (X> + D T ) < A max (£> + D T ) < 5+(H) < (2) 
Proof. Straightforward. Here the infinity norm ||-ff||oo is defined as 

\\H\\ 0O =m^\\H(vj)\\ 2 

□ 

Note that we can replace ||-H"||oo m © by an upper bound such as the one given in fl4j |. 
3. MATRIX NEARNESS CONSIDERATIONS 

Theorem 3.1. Let A = A H be any Hermitian matrix with eigendecomposition A — UAU H , with U a 
unitary and A a real diagonal matrix. Then the positive semidefinite Hermitian matrix nearest to A, both 
with respect to the Frobenius and spectral norms, is given by A + — [/max(0, A)U H . 

Proof. First we give the proof for the Frobenius norm. We need to find 

min \\X - A\\ F 

Putting X = UYU H , and exploiting the unitary invariance of the Frobenius norm, we obtain 

\\X A\\% = \\Y-A\\% = J2 + \Yu - Ad 2 

It is clear that the minimum occurs when = for i ^ j, in other words when Y is diagonal. Hence we 
obtain 

\\X - Af F = \\Y - A\\% = \Ya ~ A«| 2 

i 

It is easy to see that we must take Yu = max(0, An) and this completes the proof for the Frobenius norm. 
Note that 

mm\\X-A\\ F = I HA? 

yA,(A)<0 

For the spectral norm, it is known 0, [l(| that 



In other words, 



Now 



min \\X - A\\ 2 = infjr > : A + rl > 0} 



min \\X - Ah = max(0, -A mi r , (A)) 
x>o 



\A + -A\\ 2 - max \\,{A)\ 

Ai(A)<0 



which is zero when there are no negative eigenvalues, and — X m \n{A) when there are negative eigenvalues. □ 

Remark 3.1. From Theorem \3.1\ it is possible to find the point-wise nearest positive semidefinite matrix for 
the Hermitian matrix R{uS) = H(iu) + H(iuj) H . Obviously, if we decompose R(uj) as 

R(oj) = U(lo)A(lo)U(uj) h 

then the point-wise nearest positive semidefinite matrix is 

R+(w) = C/(w)max(A(w),0)[/(w) H 

Unfortunately, in general, the entries o/i?+(w) will not consist of rational functions and therefore cannot 
represent the transfer function of an LTI model on the imaginary axis. This problem, which in fact amounts 
to a rational approximation problem, will be addressed in the sequel. 
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4. RATIONAL APPROXIMATIONS 

Theorem 4.1. Let H(s) be passifiable , i.e., — oo < S-(H) < < 6+(H) < oo, and let R(u>) = H(iuj) + 
H(iu>) H . Let further f '(x) be a real- rational functio^ satisfying 

a > f(x) - max(x, Q) > Vx e [8-(H), 8 + (H)} (3) 

for some finite positive a. Then f(R(uj)) is positive semidefinite for all w£l Furthermore we have 

\\f(R(oj)) - -R+MHa < a VweM 

Proof. We have 

f(R(u)) - R+(u) = U{uo) {/(A(w)) - nua(A(w), 0)} f/(^) H > 

Since is positive semidefinite, the same holds for f(R(uS)). Now, since the spectral norm is unitarily 

invariant, we have 

\\f(R(uj)) - R+(oj)\\ 2 < max|/(A l H)-max(A i H,0)| 

< maxmax |/(Ai(w)) — max(Ai(w), 0)| 

max { fix) — max(x, 0)1 < a 

where the last inequality follows from the fact that all are inside the interval [S-(H), 8+(H)} . This 

completes the proof. □ 

Theorem 14 . 1 1 shows that the matrix R+(ui) can be approximated from above by the matrix f(R(w)). The 
problem is to find a suitable real-rational function f{x). We have the following : 



Theorem 4.2. Let ( n (x) = x(l + x) n /((l + x) n - 1). Then 

->Cn(x)-max(x,0)>0 Vx > -1, n = 1,2,3,... 
n 

Proof. First we prove that ( n (x) — % > for x > 0. We have 

Cn(x)-X = ^ J 

which is a positive and decreasing function for x > 0. Next we prove that Cn{x) is increasing for all x > — 1. 
This is equivalent to proving that £ n (t — 1) = — t n )/(t n — 1) is increasing for all t > 0. This is clearly 

the case for n = 1. Taking derivatives, we have 

J j-n-l 

-C„(t-l)=[n-(n + l)t + r+ 1 ] j—^ 

Now — (n + l)t + t n+1 is n when t = and 00 when t = 00. Since the derivative of n — (n + l)t + t n+1 
is (n + l)(t n — 1), the function n — (n + l)t + t n+1 attains its unique minimum (with value zero) at t = 1. 
Hence C, n {x) is increasing for all x > — 1. We therefore conclude that ( n (x) ~ max(cc,0) increases from to 
l/n in the interval [—1,0], and decreases from 1/n to in the interval [0, 00], which completes the proof. □ 

Corollary 4.1. Let H(s) be passifiable. Then the real-rational function f(x) = vC^nixjv) with v = \8-(H)\ 
satisfies the premises of Theorem \4-l\ with a = v/n. 



2 A real-rational function f(x) is a rational function assuming only real values for all real x. 
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Proof. Straightforward. 



□ 



Also, we need to find ways and means to define the matrix f(R(u:)) = f(H(iuj) + H(iu>) H ) in the whole 
s— plane and then to extract a Hurwitz stable transfer function from it. By analytical continuation, we find 
the transfer function V(s) = f (H (s) + H (— s) T ) in the entire s— plane. Since f{x) is real-rational, the transfer 
function V(s) represents the realization of a per-symmetric LTI model, i.e., satisfying V(s) = V(—s) T . This 
implies that the poles of V(s) admit the imaginary axis as symmetry axis. The following proposition 
indicates how, starting from a per-symmetric LTI model V(s) we can find a Hurwitz stable transfer function 
by additive decomposition [III EH ■ 

Proposition 4.1. Let V(s) be per-symmetric, i.e., V(s) = V{— s) T , such that V(s) has no poles on the 
imaginary axis. Then V(s) can be decomposed as V(s) = X(s) + X(—s) T , where X(s) is Hurwitz stable. 

Proof. Putting V(s) — Vq(s) + D, where Vq(s) is strictly proper and D = D T — V(oo), we can decompose 
Vo(s) uniquely into its stable and anti-stable parts, i.e., 

V (s) = X st ab(s) + X anti (s) 

Since Vo(s) is per-symmetric we have 

X s tab{s) + X ant i(s) = X sta b(~ s) T + X ant i( — s) T 

and hence X ant i(s) — X sta b{~ s) T . It follows that V(s) can be decomposed as V(s) = X(s) + X(— s) T , with 
X(s) = X s tab(s) + h D + E, where E is an arbitrary skew-symmetric matrix. It should be noted that the 
procedure is unique when the skew-symmetric matrix E is known a priori. □ 



Remark 4.1. Proposition \4-l\ assumes that V(s), in our case V(s) — f(H(s) + H(—s) T ), does not admit 
poles on the imaginary axis. By the inequality constraints 0) we know that 

a > /(Ai(w)) - nua(Ai(w), 0) > (4) 

for all eigenvalues Xi(u)) of R(uj). Since H(s) is assumed Hurwitz stable, R(ui) = H(iui) + H(iui) H cannot 
admit real poles, and hence, by the inequalities the functions /(Aj(u;)) are bounded. It follows that all 
entries ofV(ioj) = f(R(uj)) are bounded, which implies that V(s) cannot have poles on the imaginary axis. 

In the sequel we will use the Matlab® Robust Control Toolbox [l5[ routine stabproj based on the 



stable, anti-stable decomposition algorithm [12 1. 



5. TWO ALGORITHMS 

By Theorem l4.1l and Corollary |4.1l we need to find an LTI model with transfer function <p n (H(s) + H T (— s)) 
where the real-rational function 4> n {x) of denominator degree n and numerator degree n + 1 is 

/ / \ j- t i \ x{l + x/v) n 
(p n [x) = vC, n {x/v) - 



(1 + x/v) n - 1 

where v = \8—(H)\. Now it is easy to show that the following recurrence relationship holds : 

^ = 2Mx) - x - =1 ' 2 ' 3 '--- 

with 4>\(x) = x + v. 

A first algorithm ( Algorithm 1) that comes readily to the mind with Z(s) = H(s) + H T (—s) is : 
Initial value : 

Z (s) = Z(s) + vlp 
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Loop : 

for k = l to m : Z k {a) = Z k - 1 {a)(2Z k - 1 (a)-Z(a))- 1 Z k - 1 {a) 

It is seen that the associated a k upper bound at each step Z k (s), k — 0,1,..., Hi, is a k = v/2 k , and 
all Zk(iu>) are, by construction, positive semidefinite. Since the Z k {s) are all per-symmetric, we can use 
Proposition 14. 1 1 to decompose all (or only the nith one) Z k (s) in their stable and anti-stable parts as 

Z k {s) = Z s k ta \ S ) + Zl tab (-s) T 

As a last, but necessary step, we must add the skew-symmetric matrix \{D — D T ), since this matrix 
gets deleted when making the sum Z(s) = H(s) + H T (—s). In other words, the passive Hurwitz stable 
approximant G k (s) is 

G k { S ) = Zt\ S ) + \{D-D T ) 
As a very simple illustrative example take k — 0. Since Zq(s) = H(s) + H T (—s) + i/I p , we obtain easily that 

GoOO = Zt b {s) + \{D-D T ) 

(H(s) -D) + \{D + D T + ul p ) + \{D- D T ) 



which is passive by construction. In practice, Algorithm 1 has the drawback that the transfer functions 
Z k (s) in the algorithmic loop may not be minimal realizations, and hence it could happen that the stable 
anti-stable projection by means of the routine stabproj might not perform well. 

Before proposing a second algorithm, and in order to address the computational complexity of the passivated 
transfer function G(s), we want to estimate the number of poles of G(s). We suppose that f{x) is an 
irreducible real-rational function with denominator degree M and numerator degree M + 1. In this paper 
this is always the case, see also Section [51 Hence, if we suppose that all the poles are simple, we can 
decompose f(x) into partial fractions as 

M 

f(x) = a Q + /3 x + V — 

Now if the original Hurwitz stable transfer function H(s) has N poles, then the transfer function Z(s) = 
H{s) + H{—s) T has 2N poles. Also, f(Z(s)) can be written as 

M 



f{Z{s)) = a I p + (3 Z(s) + c*k(Z(s) - (3 k I P Y 



k=l 



Hence, the set of poles of f(Z(s)) is at most the union of the sets of poles of Z(s) and (Z(s) — /3fcJ p ) _1 . It 
is well known 16], that when a transfer function H(s) is such that H(oo) is invertible, then i?(s) _1 exists 
and has the same number of poles as H(s). Therefore, the number of poles of f(Z(s)), not considering 
potential cancellations, is 2N(M + 1). Finally, after the stable anti-stable decomposition, this number is 
to be divided by two, to yield N(M + 1) poles for the final passivated transfer function G(s). Of course 
the number N(M + 1) is only an estimate, since pole-zero cancellations can occur. If for some reason, the 
number of poles of the explicitly proved passive transfer function G(s) appears to be unacceptable high, a 



final judiciously chosen passivity preserving model order reduction step [17H2 1| can be applied. 

Hence, in order to find a workable algorithm, we have to find the partial fractions decomposition of f(x) 

^>n( x ) = v Cn(x/v). If we restrict ourselves to even n = 2m > 2, we have the partial fraction expansion 

( 2m (x) =x + - —— + Ke V 

777, I *n -\- J. < ^ 



2 x + 1 - e vik l m 

k—1 
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Hence, with f(x) = vC,im(%iv) we have 



rn — l 



f(Z(s)) = Z(s) + — [Z(s) + 2vl p \- x + 5Re ^ L?™ h l m - e *ik/m\ [z(s) + ^ _ ^ik/m^-i (5) 



k=i 



Algorithm 2 performs the state space addition ([5]) as is, i.e., we add the realizations of Z(s), [y 2 / m)[Z (s) 
2vl.p}~ 1 , etc., to obtain f(Z(s)). The explicit state space form for the terms 



3?e 



in formula ([5]) is obtained by the state space technique described in the Appendix. Finally, the stable 
anti-stable projection yields the passivated transfer function G(s). 

5.1. Numerical Examples 

We will consider only reciprocal non-passive systems, i.e., systems with H(s) — H(s) T , as these systems 
are representative of LTI systems satisfying the electromagnetic condition known as Lorentz reciprocity [22[ . 
Of course the theory also remains valid for non-reciprocal LTI systems. Since for reciprocal systems R(cj) 
is real and even, this explains why the plots in the sequel only show values for non-negative frequencies. 

5.1.1. First example 

As a first example we take the SISO Hurwitz stable non-passive transfer function 



H(s) 



7.2s 4 + 47-Ol.s 3 + 230.8s 2 + 536.6s + 587.1 
3.2s 4 + 32.61s 3 + 43.63s 2 + 117.5s + 104.3 



(6) 



We use the approach of Algorithm 1 with n\ = 2. The passivated approximation G(s) has a non-minimal 
realization with 65 poles which are reduced to 20 by the routine minreal [l6j]. The real and imaginary parts 
of the original transfer function H(s) vs. the passivated transfer function G(s) are shown in Figs [T] and [2l 



Real part of the transfer function 




Original 
- Passivated (Alg. 1, n 1= 2) 



Figure 1: Real part of G(s) vs. H(s) 



5.1.2. Second example 

As a second example we take the SISO Hurwitz stable minimum phase non-passive transfer function 

= (s + l)(s + 3)(s + 90)(s + 95)(s + 100) 

U (s + 25)(s + 35)(s + 38)(s + 180)(s + 185) 1 ' 
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Imaginary part of the transfer function 




CO 



Figure 2: Imaginary part of G(s) vs. H(s) 



We use the approach of Algorithm 2 with m = 5. The passivated approximation G(s) has a realization with 
50 poles. The real and imaginary parts of the original transfer function H(s) vs. the passivated transfer 
function G(s) are shown in Figs [3] and SJ 

Real part of the transfer function 
0.16 I . 1 , 1 1 , 1 1 




-0.02- 

-0.04 ' ' ^— - — 1 1 1 1 

5 10 15 20 25 30 35 40 

CO 



Figure 3: Real part of G(s) vs. H(s) 



5.1.3. Third example 

As a third example we take the 2x2 MIMO Hurwitz stable non-passive transfer function 

o ,12 2s + 10 



H(s) 



s z + 3s 4 
2s + 10 



3s + 2 



(8) 



We use the approach of Algorithm 2 with m = 4. The passivated approximation G(s) has a realization 
with 48 poles. Fig. [S] plots the values of X m i n {G{iuj) + G{iuj) H ) vs. \ min (H(iuj) + H(iuj) H ). To show 
the nearness of the original and passivated transfer functions H(s) and G(s), we plot the relative error 
||G(iw) - H(w)\\ 2 /\\H(iL>)\\2 in Fig. 1 



Imaginary part of the transfer function 
0.35 I 1 1 1 1 , r- 




Figure 4: Imaginary part of G(s) vs. H(s) 



Minimum eigenvalue of R(to) 




CO 



Figure 5: Minimum eigenvalues of passivated vs. original transfer functions 



6. MINIMAX ALGORITHM 

The starting point for finding a passive approximant is to find a real-rational function f(x) that satisfies 

a > f(x) - max(i, 0) > Vx G [-a, b] (9) 

where a = —6-(H) = \8-(H)\ and b — 8+(H). Since max(cc,0) = (|x| + x)/2, this can be written as 

a > 2f(x) - x - a - \x\ > -a Va; G [-a, b] (10) 

Putting r(x) = 2/(x) — x — a, and since our aim is to find the smallest positive a such that (If 0[) is satisfied, 
it is seen that we must find the rational minimax or Chebyshev approximant, i.e., 

min max \r(x) — \x\\ 

r x£[— a.b] 



Let us first treat the case a = b = 1, which is well-documented in the literature 23M25j. Since \x\ is even and 
the interval [—1, f] is symmetric with respect to 0, it is clear that r(x) must be an even rational function, i.e., 
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||G(i to) - H(i to) || 2 /||H(i to)|| 2 




to 



Figure 6: Relative error between passivated and original transfer functions 



r(x) — p{x 2 ). If we take p(t) irreducible with numerator and denominator of exact degree n, the minimax 
problem can be reformulated as: 

min max \p(t) — (11) 

p 0<t<l 

Calling E n the value obtained by the minimax problem (|11[) , it is clear that at the minimum we must have 

E n > p(t) -Vt> -E n for < t < 1 (12) 

Furthermore, the Remes condition (2f| [26[ requires that there are exactly 2n + 2 point tk inside [0, 1] where 
the equality 

Vh-p{t k ) = (-l) k E n fc = l,2,...,2n + 2 

is satisfied. This allows an iterative approach [25[ to find the optimal E n and p(t). The poles and zeros 
of p(t) are all simple and intertwined on the negative real axis |27| . It follows that in general p(t) can be 
written as 

n 

p(t) = a - V — ^~ 

k—l 

where all a,k-,bk are positive. For n — 4 the coefficients ak,bk with &q = E n are given in Table [TJ Fig. [7] 



Table 1: Coefficients for the function p(t) for n = 4 








1 


2 


3 


4 




2.6397296257 


1.4034219887xl0- B 


0.0003730797 


0.0290141901 


5.6266532592 




0.0007365636 


0.0000917473 


0.0049831021 


0.1014048457 


2.4866930733 



shows the approximation error p(t) — \ft and the equioscillation property. Note that the asymptotic formula 
of E n is known [28], i.e., we have E n w 8e _7rv ^™ for n — > ex). 
Formula (1X21) implies 

2£„ > p(x 2 ) + £„ - |a;| > for - 1 < x < 1 

or 

S„ > P(a;2) y + g " - max(0,s) > for - 1 < x < 1 (13) 
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t 



Figure 7: Minimax approximation error for n = 4. 



For a = b = 1, the best rational function f(x) satisfying © is therefore f(x) = ^(p(x 2 ) + x + E n ) with 
a = E n . It should be noted that f(x) has numerator degree 2n + 1 and denominator degree 2n. The case of 
the general interval [— a, b] instead of [—1,1] is treated by the following 

Theorem 6.1. Let a, b > and f{x) a real-rational function such that 

a > f(x) — max(i, 0) > for — 1 < x < 1 

Then the real-rational function 



fa,b(x) 



x(b — a) + lab 



f 



x(b + a) 
x(b — a) + 2ab 



is such that 



a max(a, b) > f a ,b(x) — max(x, 0) > for — a < x < b 



Proof. The bilinear transformation g{x) — x{b + a)/(x(b — a) + 2ab) maps the interval [—a, b] onto the 
interval [—1, 1]. Moreover, the linear function x(b — a) + 2ab is positive over [—a, b] since it is positive at the 
endpoints. Hence 

a > f(g( x )) ~ max (g( x ): 0) > for — a < x < b 



implying 



x(b — a) + 2ab 

a ; > 

a + b 



x(b — a) + 2ab 



f(g(x)) — max(a;, 0) > for — a < x < b 



which completes the proof. Note that, if the denominator degree of f(x) is m and the numerator degree is 
m + 1, then the same holds for f a ,b(x). □ 

In light of formula (flU]) . we take f(x) = ^(p(x 2 ) + x + E n ) and a = E n . The function f a ,b(x) can be 
conveniently written as 

fa,b( X ) = ( TX + K )f 



where 



6 + (H)- \8-(H)\ 
5+(H) + \5-(H)\> 
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, 6 + (H)\S-(H)\ 
'S + (H) + \S.(H)\ 



For the function f l (x) — l/(x 2 + £), the partial fraction expansion of /r b (x) is given by 

(rx + nf t 3 t 2 (3 + t 2 £) m f »/(t,k,/) 

— •■ 1 Me 



x 2 +£(tx + k) 2 " 1 + t 2 £ (1 + t 2 £) 2 



-t(T,K,t) 



where 



Kyf£ 



-VI 



r)(r, k,£) 



£(t,M) £ 
k£ 2 



Hence for the function fix) = ■k(p(x 2 ) + x + E n ), the transformed function f a ,b(x) can be written as 



fa,b( X ) 



(E n + o, )(tx + k) a kfa%(z 



k=l 



(14) 



The partial fraction expansion of (fTl)) is the key of Algorithm 3, since we ultimately have to calculate 
fa,b(Z(s)), where Z(s) — H(s) + H(—s) T . The linear terms of (|14p all add up to the compound linear term 



n , 

x + (E n + Oq){tx + k) - ^2 a k I x- 

fc=i ^ 



r 3 r 2 (3 + T 2 b k ) 

+ T 2 b k +K (1+T 2 b k ) 2 



leading to a linear term fab ear (Z(s)) — k\Z(s) + k2l p - The remaining terms, obtained by evaluating 

Me {^t, k, b k )(Z{s) - £(t, «, fofc)^)- 1 } 

are obtained by the state space technique described in the Appendix. Finally, as in Algorithm 2, the stable 
anti-stable projection of f a ^(Z(s)) is performed in order to obtain the passivated transfer function G(s). 

6.1. Numerical Examples 
6.1.1. First example 

As our first example we again take the SISO Hurwitz stable minimum phase non-passive transfer function 
([7]) , but here we use Algorithm 3 with n = 4 and the coefficients of Table Q] The passivated approximation 
G(s) has a realization with 45 poles. The real and imaginary parts of the original transfer function H(s) 
vs. the passivated transfer function G(s) are shown in Figs [5] and [HI It is seen by comparing with Figs [3] 
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Figure 8: Real part of G(s) vs. H(s) 



and [4] that the approximation is better, while requiring 5 poles less. 
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Figure 9: Imaginary part of G(s) vs. H(s) 



6.1.2. Second example 

For the second example we again take the MIMO Hurwitz stable non-passive transfer function (O, but here 
we use Algorithm 3 with n = 4 and the coefficients of Table [T] The passivated approximation G(s) has a 
realization with 46 poles. Fig. rXU]plots the values of X m i n (G(iuj) + G(iuj) H ) vs. X m i n (H(iuj) + H(iuj) H ). To 
show the nearness of the original and passivated transfer functions H(s) and G(s), we plot the relative error 
\\G(iuj) — H(iuj)\\2/\\H(iuj)\\2 in Fig. [TT] It is seen by comparing with Figs [5] and H2 that the approximation 
is more or less similar, but requires 2 poles less. 
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Figure 10: Minimum eigenvalues of passivated vs. original transfer functions 



7. CONCLUSION 

We have presented a new global passification approach towards finding a passive transfer function G(s) 
that is nearest in some well-defined matrix norm sense to a given non-passive transfer function H(s). It 
is shown that the key point in constructing the nearest passivated transfer function G(s), is to find a 
good rational approximation to the well-known ramp function over an interval defined by the minimum and 
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Figure 11: Relative error between passivated and original transfer functions 

maximum dissipation of the given non-passive transfer function H(s). It is also shown that in the Chebyshev 
or minimax sense this requires finding a rational Chebyshev approximation of the square root function over 
the unit interval. The proposed algorithms rely strongly on the stable anti-stable projection of a given 
transfer function. Five pertinent examples, both SISO and MIMO, are given to show the accuracy and 
relevance of the proposed algorithms. 

8. APPENDIX 

Suppose we have the real-rational transfer function H(s) = C(sl n — A)^ 1 B + D, and we need to evaluate 
the transfer functioijf] 

Zleri{H(8) - SI,)- 1 

where 77 and £ are complex numbers. Suppose also that D — £I p is invertible. Putting = D — £I p , we 
have [l6| that the complex state space transfer function rj(H(s) — is given by C(sl — A)^ 1 !} + D 

where 

A = A — BD^C, B = qBD-T 1 , C = -D^C, D = tjD^ 1 
In state space form we have 

x = Ax + Bu 
y = Cx + Du 

The input u is real, but the output y is complex. Putting y = yi + iy-2, it is clear that we are only interested 
in yi as output. Decomposing all complex vectors and matrices in their real and imaginary components, we 
obtain 

xi+ix 2 = (Ai + iA 2 )(xi + ix 2 ) + (Bi + iB 2 )u 
IJ1+W2 = (Ci + iC 2 )(xi + ix 2 ) + (Di + iD 2 )u 

Hence the state space equations with u as input and y\ as output are simply : 

x\ = A\X\ — A 2 x 2 + B\u 
x 2 = A 2 xi + Aix 2 + B 2 u 
yi = C\x\ - C 2 x 2 + Diu 



||G(i 03) - H(i to) || 2 /||H(i co)|| 2 




CO 



Considering s = d/dt as a real operator. 



14 



References 



[1] B. Anderson and S. Vongpanitlerd, Network Analysis and Synthesis, NJ: Prentice-Hall, 1973. 

[2] S. Grivet-Talocia, Passivity enforcement via perturbation of Hamiltonian matrices, IEEE Trans. Circuits Syst. I 51 (9) 

(2004) 1755-1769. 

[3] B. R. Andrievskii, A. L. Fradkov, The passification method in problems of adaptive control, estimation, and synchroniza- 
tion, Autom. Remote Control 67 (11) (2006) 1699-1731. 

[4] C. J. Damaren, H. J. Marquez, A. G. Buckley, Optimal strictly positive real approximations for stable transfer functions, 
IEE Control Theory and Applications 143 (6) (1996) 537-542. 

[5] B. Dumitrescu, Parameterization of positive-real transfer functions with fixed poles, IEEE Trans. Circuits Syst. I 49 (4) 
(2002) 523-526. 

[6] C. P. Coelho, J. Phillips, L. M. Silveira, A convex programming approach for generating guaranteed passive approximations 
to tabulated frequency-data, IEEE Trans. Computer- Aided Design 23 (2) (2004) 293-301. 

[7] D. Saraswat, R. Achar, M. S. Nakhla, Global passivity enforcement algorithm for macromodcls of interconnect subnetworks 
characterized by tabulated data, IEEE Trans. VLSI Systems 13 (7) (2006) 819-832. 

[8] Y. Tanji, H. Kubota, Passive approximation of tabulated frequency data by Fourier expansion method, in: Proc. ISCAS, 
vol. 6, 2005, pp. 5762-5765. 

[9] N. J. Higham, Matrix nearness problems and applications, in: Applications of Matrix Theory, University Press, 1989, pp. 
1-27. 

[10] P. R. Halmos, Positive approximants of operators, Indiana Univ. Math. J. 21 (1971/72) 951-960. 

[11] B. Kagstrom, P. VanDooren, A generalized state-space approach for the additive decomposition of a transfer matrix, J. 

Numcr. Linear Algebra Appl. 1 (2) (1992) 165-181. 
[12] M. G. Safonov, E. A. Jonckheere, M. Vermaj, D. J. N. Limebeer, Synthesis of positive real multivariablc feedback systems, 

International Journal of Control 45 (3) (1987) 817-842. 
[13] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan, Linear Matrix Inequalities in System and Control Theory, SIAM 

Studies in Applied Mathematics 15, Philadelphia, PA, 1994. 
[14] S. Boyd, V. Balakrishnan, P. Kabamba, A bisection method for computing the norm of a transfer matrix and related 

problems, Math. Control Signals Systems 2 (3) (1989) 207-219. 
[15] G. Balas, R. Chiang, A. Packard and M. Safonov, Robust Control Toolbox User's Guide, Version 3, The MathWorks, Inc, 

2005. 

[16] Control System Toolbox User's Guide, Version 4, The MathWorks, Inc, 1998. 

[17] X. Chen, J. T. Wen, Positive realness preserving model reduction with norm error bounds, IEEE Trans. Circuits 

Syst. I 42 (1) (1995) 23-29. 
[18] L. Knockacrt, A note on strict passivity, Systems Control Lett. 54 (9) (2005) 865-869. 

[19] D. C. Sorensen, Passivity preserving model reduction via interpolation of spectral zeros, Systems Control Lett. 54 (4) 

(2005) 347-360. 

[20] A. C. Antoulas, A new result on passivity preserving model reduction, Systems Control Lett. 54 (4) (2005) 361-374. 
[21] L. Knockacrt, T. Dhaene, F. Ferranti, D. DeZutter, Model order reduction with preservation of passivity, non-expansivity 

and Markov moments, Systems Control Lett. 60 (1) (2011) 53—61. 
[22] L. Knockaert, D. DeZutter, On the complex symmetry of the Poincare-Stcklov operator, Progress in Electromagnetic 

Research B. 7 (2008) 145-157. 
[23] D. J. Newman, Rational approximation to |x|, Michigan Math. J. 11 (1964) 11-14. 
[24] L. Brutman, E. Passow, On rational approximation to |x|, Constr. Approx. 13 (1997) 381-391. 

[25] R. S. Varga, A. Ruttan, A. D. Karpenter, Numerical results on the best uniform rational approximations of the function 

|x| on the interval [-1, 1], Mat. Sb. 182 (11) (1991) 1523-1541. 
[26] A. Ralston, Rational Chebyshev approximation by Rcmes algorithms, Numer. Math. 7 (1965) 322-330. 
[27] H. P. Blatt, A. Iscrlcs, E. B. Saff, Remarks on the behaviour of zeros of best approximating polynomials and rational 

functions, in: Algorithms for approximation, Inst. Math. Appl. Conf. Ser. New Ser., 10, Oxford Univ. Press, 1987, pp. 

437-445. 

[28] H. R. Stahl, Best uniform rational approximation of x a on [0, 1], Acta Math. 190 (2003) 241-306. 



15 



